clear
use "$maindir/Results/Model_Outputs/ihwap_wedge.dta"

*** covariate means with state sample
* sample to be used for descriptive statistics
reg tau_ml i.treated if twoyears_prepost==1 & meterdays_monthly!=. & regsample==1, nocons
gen descrsample = e(sample)

count if ww_treat==1 & ProgramYear>=2009 & descrsample==1
sca homes_util = `r(N)'

replace sqfeet = . if sqfeet==0
winsor2 sqfeet, cuts(0.5 99.5) trim replace

gen pct_BlowerReduc = BlowerReduc/Blower_Pre

gen atticR = Attic_RValue
replace atticR = . if atticR==0

gen liheap = 0 if liheapemergency=="False"
replace liheap = 1 if liheapemergency=="True"

gen tabbuilddateX = 0 if builddate!=.
replace tabbuilddateX = 1 if tab_builddate11==1 | tab_builddate12==1 | tab_builddate13==1

gen decade1 = tab_builddate1==1
gen decade2 = tab_builddate2==1 | tab_builddate3==1 | tab_builddate4==1
gen decade3 = tab_builddate5==1 | tab_builddate6==1 | tab_builddate7==1
gen decade4 = tab_builddate8==1 | tab_builddate9==1 | tab_builddate10==1
gen decade5 = tabbuilddateX==1


foreach x in female renter white black hispanic nativeamerican otherrace ///
		haselderly hasminor liheap pct_BlowerReduc ///
		nstories decade1 decade2 decade3 decade4 decade5 ///
		{
		replace `x' = `x'*100
}


matrix input ttest_mat = (.,.,.,.)
foreach x in Real_income1000 occupants Age female renter white black hispanic nativeamerican otherrace ///
		haselderly hasminor liheap Blower_Pre Blower_Post BlowerReduc pct_BlowerReduc ///
		atticR sqfeet bedrooms nstories decade1 decade2 decade3 decade4 decade5 ///
		Real_total_cost Real_tot_actAirCon Real_tot_actAirSeal Real_tot_actAttic  Real_tot_actBaseload Real_tot_actDoor ///
		Real_tot_actFoundation Real_tot_actFurnace Real_tot_actGeneral Real_tot_actHealSfty ///
		Real_tot_actWallIns Real_tot_actWtHtr Real_tot_actWindow {

sum `x' if ww_treat==1 & ProgramYear>=2009 & descrsample==1 , detail
sca mean`x' = `r(mean)'
sca sd`x' = `r(sd)'
sca min`x' = `r(min)'
sca max`x' = `r(max)'
matrix mat`x' = (mean`x' , sd`x', min`x', max`x')
matrix ttest_mat = (ttest_mat \ mat`x')
}
matrix homes = (homes_util,.,.,.)
matrix ttest_mat = (ttest_mat \ homes)
matrix ttest_mat = (ttest_mat \ .,.,.,.)

foreach x in gas_mmbtu electric_mmbtu tmin tmax precip {
	sum `x' if ProgramYear>=2009 & descrsample==1 , detail
	sca mean`x' = `r(mean)'
	sca sd`x' = `r(sd)'
	sca min`x' = `r(min)'
	sca max`x' = `r(max)'
	matrix mat`x' = (mean`x' , sd`x', min`x', max`x')
	matrix ttest_mat = (ttest_mat \ mat`x')
}
count if regsamp==1 & ProgramYear>=2009 
sca obs_month = `r(N)'
matrix obs = (obs_month,.,.,.)
matrix ttest_mat = (ttest_mat \ obs)

matrix colnames ttest_mat = "Average" "Standard Deviation" "Min" "Max" 

matrix rownames ttest_mat = "" "Income ($/1000)" "N Occupants" "Householder Age" "Female Householder (%)" ///
		"Renter (%)" "White (%)" "Black (%)" "Hispanic (%)" "Native American (%)" "Other Race (%)" ///
		"Seniors 65+ (%)" "Children Under 18 (%)" "LIHEAP (%)" "Blower Door Pre (CFM50)" ///
		"Blower Door Post (CFM50)" "Blower Door Reduced (CFM50)" "Percent Blower Door Reduced (%)" ///
		"Attic R-Value" "Floor Area (sqft)" "N Bedrooms" "Has Multiple Stories (%)" "Built Pre-1900 (%)" ///
		"Built 1900-1929 (%)" "Built 1930-1959 (%)" "Built 1960-1989 (%)" ///
		"Built 1990-Present (%)" "Total" "Air Conditioning" "Air Sealing" ///
		"Attic" "Baseload" "Door" "Foundation" "Furnace" "General" "Health and Safety" "Wall Insulation" ///
		"Water Heater" "Window" "Number of Homes" "" "Monthly Gas Usage (MMBtu)" "Monthly Elec. Usage (MMBtu)" ///
		"Min. Temperature (C)" "Max. Temperature (C)" "Precipitation (in)" "Number of Observations"

esttab matrix(ttest_mat, ///
	fmt("%9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.0fc %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.0fc")) ///
	using "$maindir/Results/Tables/descriptives.tex", replace
	
	

************************************************	
**** GRAPH ILLUSTRATING WEATHERWORKS BIAS
sort Household meterend

** identify homes with a full year of pre or post treatment data
by Household : gen days_diff = meterend - meterend[_n-12]
by Household : gen days_diff_post = meterend[_n+12] - meterend
by Household : gen year_home_pre = 1 if days_diff>=360 & days_diff<=370 & treated==0 & treated[_n+1]==.
by Household : egen full_year_pre = max(year_home_pre)
drop year_home_pre
by Household : gen year_home_post = 1 if days_diff_post>=360 & days_diff_post<=370 & treated==1 & treated[_n-1]==.
by Household : egen full_year_post = max(year_home_post)
drop year_home_post

** total energy consumption for a full year pre or post treatment
gen mmbtu_year_pre = total_mmbtu + total_mmbtu[_n-1] + total_mmbtu[_n-2] + total_mmbtu[_n-3] + total_mmbtu[_n-4] + total_mmbtu[_n-5] + ///
	total_mmbtu[_n-6] + total_mmbtu[_n-7] + total_mmbtu[_n-8] + total_mmbtu[_n-9] + total_mmbtu[_n-10] + total_mmbtu[_n-11] ///
	if treated==0 & treated[_n+1]==. & full_year_pre==1

gen mmbtu_year_post = total_mmbtu + total_mmbtu[_n+1] + total_mmbtu[_n+2] + total_mmbtu[_n+3] + total_mmbtu[_n+4] + total_mmbtu[_n+5] + ///
	total_mmbtu[_n+6] + total_mmbtu[_n+7] + total_mmbtu[_n+8] + total_mmbtu[_n+9] + total_mmbtu[_n+10] + total_mmbtu[_n+11] ///
	if treated==1 & treated[_n-1]==. & full_year_post==1

gen usageratio_pre = (ExistingConsumption/1000000)/(mmbtu_year_pre)
gen usageratio_post = (ProjectedConsumption/1000000)/(mmbtu_year_post)

foreach x in usageratio_pre usageratio_post {
winsor2 `x', replace cuts(0.5 99.5)
}

****** Plotting the graphs of ratio Wx Usage/True Usage
**** True usage on X-axis
gen tot_year_prebins = .
forval i = 40(20)180 {
replace tot_year_prebins = `i' if mmbtu_year_pre>=`i' & mmbtu_year_pre<`i'+20
}
replace tot_year_prebins = 20 if mmbtu_year_pre<40
replace tot_year_prebins = 200 if mmbtu_year_pre>=200 & mmbtu_year_pre!=.
bysort Household : egen total_year_prebins = max(tot_year_prebins)
drop tot_year_prebins

reg usageratio_pre ibn.total_year_prebins, nocons
est sto usageratio_pre
reg usageratio_post ibn.total_year_prebins, nocons
est sto usageratio_post

coefplot (usageratio_pre, mcolor(gs4) msymbol(Th) ciopts(lcolor(gs4))) (usageratio_post, mcolor(black) ciopts(lcolor(black))), ///
		offset(0) vertical xlabel(1 "<40" 2 "40-60" 3 "60-80" 4 "80-100" 5 "100-120" ///
		6 "120-140" 7 "140-160" 8 "160-180" 9 "180-200" 10 ">=200", labsize(small) angle(45)) ///
		ylabel(0(1)6, glcolor(gray%10)) ///
		xtitle("Average Yearly Pre-Treatment Energy Consumption (MMBtu)") ///
		ytitle("Ratio = Modeled/Actual Energy Usage") legend(order(2 "Pre-Treatment" 4 "Post-Treatment")) ///
		rename(*.total_year_prebins = .total_year_postbins) ///
		yline(1, lcolor(black)) graphregion(color(white)) bgcolor(white)
graph export "$maindir/Results/Figures/ratio_totalbins.pdf", replace
